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PREFACE 

An  adaptive  controller  is  one  which  has  the  capability  of 
learning  about  unknown  aspects  of  a  system  being  controlled  and  then 
modifying  its  control  regime  in  an  effort  to  improve  the  quality  of 
the  control  exerted. 

This  paper  is  devoted  to  a  mathematical  formulation  and  com¬ 
putational  solution  of  the  problems  of  system  identification  and 
the  determination  of  unmeasurable  state  variables  on  the  basis  of 
observations  of  a  process,  two  topics  of  central  importance  in  the 
design  of  adaptive  controllers. 

The  approach  suggested— based  on  the  theory  of  quasilineari¬ 
zation- -is  an  outgrowth  of  continuing  RAND  research  on  the  computa¬ 
tional  solution  of  multi -point  boundary-value  problems.  The  paper 
should  be  of  interest  to  control  engineers  and  numerical  analysts. 
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SUJWAOT 


The  dynamical  equations  for  a  system  are  prescribed  in  the 

form 

(1)  x  =  g(x,«),  x(o)  =  c, 

where  the  system  parameter  vector  a  and.  the  initial  vector  c  are  un¬ 
known.  Noisy  observations  on,  e.g.,  the  first  component  of  the  state 
vector,  x^(t),  are  made  at  times  t^, 

(2)  « 1)^,  i  =  1,2, 

where  b^  is  the  observation  at  time  t^ .  It  is  desired  to  find  an 
initial  vector  c  and  a  system  parameter  vector  or  which  minimize  the 
sum  of  the  squares  of  the  deviations 

H  2 

(3)  5  '  ^  {Xl(tl>  '  ”l}  ' 

An  effective  computational  scheme,  based  on  the  notion  of 
quasilinearization,  is  suggested.  The  utility  of  the  method  is 
illustrated  by  considering  a  process  described  by  the  Van  der  Pol 
equation,  in  which  the  system  parameter  and  the  initial  velocity- 
are  to  be  determined  on  the  basis  of  observations  of  the  displace¬ 
ments  at  various  times.  A  FORTRAN  program  is  provided. 

These  considerations  provide  a  general  approach  to  the  problems 
of  system  identification  and  the  determination  of  unmeasurable  state 
variables,  two  matters  of  prime  importance  in  the  design  of  adaptive 


controllers 
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I.  INTRODUCTION 

Much  of  the  modem  theory  of  control  processes  is  devoted  to 
the  minimization  of  functionals  of  the  form 

T 

(1)  J[y]  =  f  h(x,y)dt  , 

v 

o 

where  the  state  vector  x,  which  is  of  dimension  R,  and  the  control 
vector  y  are  related  "by  the  dynamical  equation 

(2)  x  =  g(x,y)  ,  x(0)  =  c  . 

Though  much  progress  has  been  made  in  the  analytical  and  computa¬ 
tional  treatment  of  such  problems  [1,2],  much  remains  to  be  done. 

In  particular,  it  must  be  recognized  that  frequently  it  is  not  pos¬ 
sible  to  measure  all  of  the  components  of  the  state  vector,  and  the 
dynamical  equation,  Eq.  (2),  may  not  be  known  precisely  by  the  con¬ 
troller.  Wider  these  circumstances  effective  control  is  more  diffi¬ 
cult  to  achieve,  but  not  impossible.  Some  of  the  progress  which  has 
been  made  in  treating  these  matters  is  discussed  in  [2,3].  In  addi¬ 
tion,  various  questions  arise  concerning  the  choice  of  the  index  of 
performance  in  Eq.  (l),  and  some  are  treated  in  [4]. 

In  this  Memorandum  our  aim  is  to  show  how  many  control  problems 
involving  unmeasurable  state  variables  and  unknown  system  dynamics 
lead,  from  the  mathematical  viewpoint,  to  nonlinear  multi-point  bound¬ 
ary  value  problems,  and  also  to  show  how  the  quasilinearization  tech¬ 
nique,  discussed  in  Refs.  [5,6,7],  leads  to  an  efficient  computational 
mode  of  solution.  In  effect,  a  procedure  for  the  design  of  a  broad 
class  of  adaptive  controllers  is  given. 


2 


II.  FORMULATION 

Let  us  suppose  that  the  exact  dynamical  equation  is  unknown  to 
the  controller  but  that  its  general  form  is  known,  the  only  unknowns 
being  a  vector  of  system  constants  a.  The  system  equations  are  then 
assumed  to  be 

(1)  x  =  g(x,y,or)  . 

In  addition  we  assume  that  at  certain  times  sane  of  the  components 
of  the  state  vector  are  measured  by  sensing  equipment  and  the  con¬ 
troller  is  apprised  of  these  measurements.  The  basic  problem  of  feed¬ 
back  control  is  the  determination  of  the  optimal  choice  of  the  control 
vector  y  for  any  set  of  circumstances  in  which  the  controller  may 
find  itself. 

We  shall  treat  this  problem  by  seeking  to  determine  the  values 
of  the  parameters  in  Bq.  (l)  and  the  solution  of  Eq.  (l)  which  are 
in  best  agreement  with  the  measurements,  which  reduces  the  problem 
to  the  classical  one  mentioned  in  Section  I. 

Ifore  precisely,  we  wish  to  determine  the  vector  a  and  a  complete 
set  of  initial  conditions  c  so  that  the  solution  of  the  system 

(2)  x  =  g(x,a)  ,  x(0)  =  c, 

is  in  closest  agreement  with  the  observations.  For  ease  of  descrip¬ 
tion  let  us  suppose  that  observations  of  the  first  component  of  the 
state  vector  are  available  at  times  t^  i=l,2, ...,N.  We  denote  the 
observed  value  at  time  ti  by  b^.  Our  aim  is  to  determine  the  system 
parameter  vector  a  and  the  initial  state  vector  c  that  minimize  the 
sum  of  the  squares  of  the  deviations 
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where  x1(t±)  is  the  first  component  of  the  vector  x  evaluated  at 
t  “  V  **  WBS  “entioned  earlier,  once  the  minimizing  system  vector 
and  initial  s^te  vector  are  determined,  we  shall  consider  them  to  be 
the  "actual"  system  and  initial  state  vectors  and  control  the  system 
on  that  basis.  As  the  process  continues,  the  identification  and  de- 
termination  procedures  are  to  be  repeated  from  time  to  time,  so  that 
adaptation  is  possible. 
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III.  A  PRELIMIHART  SIMPUPICATION 


It  is  inconvenient  to  consider  the  system  vector  a  and  the  initial 
state  vector  c  to  be  different  types  of  vectors.  Let  us  consider  or 
to  be  part  of  the  state  vector,  where  a  is  subject  to  the  equation 

(1)  a  =  0  . 


Then  cur  object  is  to  determine  the  initial  vectors  er(0)  and  x(0)  in 
such  a  way  that  we  minimize  the  sum 


(2) 


S  =  l  (x1(ti)  -  b^.) 
i=l 


) 


where 

fx  =  g(x,Qf) 

(3) 

Q  • 

11 

O 

But  this  is  equivalent  to  the  problem  of  minimizing  the  sum  S,  subject 
to  the  condition 

(*0  x  =  g(x),  x(0)  =  c  , 

where  the  minimization  is  over  the  initial  vectors  c,  through  a  re¬ 
interpretation  of  x,g(x), 


and  c. 
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IV.  QIXAS ILINEARIZATI ON 


Consider  the  sequence  of  vector  functions  x^(t),  x^(t),  ..., 

(2) 

xv  '(t),  ...  defined  in  this  recursive  manner.  The  vector  function 
x(°)(t)  is  chosen  on  the  interval  0  <  t  s  tN.  After  the  function 
x^(t)  is  determined,  the  vector  function  x^(t)  is  taken  to  be  the 
solution  of  the  linearized  system  of  differential  equations 


(1) 


•(k+D 

i 


)  + 


*  9gi(x(k)) 


i  -  •  •  •  jR  • 

The  initial  conditions  are  to  be  selected  so  as  to  minimize  the  sum 

(2)  Sx  -  l  (4k+1)  Ct±)  -  . 

1=1 

Then,  in  many  instances  of  practical  importance,  [5,7],  the  sequence 
of  vectors  x^(t)  vill  converge  quadratically  to  the  desired  solution 
of  the  problem  of  the  previous  section,  and  the  vectors  xv  '(0)  will 
converge  to  the  desired  minimizing  initial  vector.  Other  applications 
of  this  idea  to  orbit  determination  [8]  and  system  design  [9]  have 


been  made 
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V.  COMPUTATIONAL  CONS  EjERATIOJB  -  I 

Since  x^+1^(t)  is  the  solution  of  a  linear  system  of  differen¬ 
tial  solutions,  it  may  be  represented  in  the  form 

a)  x(k+1)(t)  =  P(t)  Cih(i)Ct)  , 

i=l 

where  the  vector  function  p(t)  is  the  particular  solution  of  Eg..  (4.1) 
subject  to  the  initial  conditions 

(2)  p(o)  =  0  , 

and  the  vector  h^(t)  is  the  solution  of  the  homogeneous  form  of 
Eq.  (4.1),  the  initial  vector  having  the  i  row  unity  and  the  others 
zero.  These  vectors,  p(t)  and  h^(t),  are  all  considered  to  he 
determined  computationally  on  the  interval  0£t£t^.  The  scalars 
c^,  i=l,2, ...,R,  are  those  that  minimize  the  sum 

(3) 

j=i  i=i 

and  may  be  determined  from  the  normal  equations  [10] 

w  7  favi  -  4iK“,(V  ■ 0  > 

j=i  i=l 

m  —  1,2, ... ,R  . 

In  this  way  the  problem  of  determining  the  optimal  initial 
vector  is  made  to  depend  upon  the  numerical  solution  of  systems  of 
linear  initial-value  problems  and  of  linear  algebraic  equations. 

Let  us  next  consider  a  special  case. 
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VI.  THE  VAN  PER  POL  EQUATION 

Suppose  that  a  process  is  described  by  Van  der  Pol's  equation 

[11,12] 

(1)  x  =  u 

•  2 

U  =  -  A.(x  -  l)u  -  X  , 

where  now  x  is  a  scalar  and  X  is  an  unknown  constant.  We  suppose 
that  the  following  three  observations  of  the  displacement  x  have 
been  made 

(2)  x(4)  =  -  1.8084-3 
x(6)  =  -  1.63385 
x(8)  =  -  1.4o456  . 


We  wish  to  determine  the  unknown  value  of  X  and  both  x  and  u  for 

t  =  4. 

We  consider  the  system  of  equations 

(3)  x  =  u 

u  =  -  x(x2-  l)u  -  x 

X  =  0 

and  wish  to  determine  x(4),  u(4)  and  X(4-)  so  that  Eqs.  (2)  will  be 
satisfied.  In  this  instance  it  is  not  necessary  to  use  the  method 
of  least  squares  since  only  three  observations  are  given. 

To  obtain  an  initial  approxi-aation  we  observe  that 

(4)  *(§)_:- w  .o87  w  x(4) 
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(5)  *(fl>  ~2^6K  .li4«,x(6) 

(6)  x(6)  -  ^  .oi4  ps  x(4.o) 

In  view  of  Eq.  (l)  ve  are  led  to  the  initial  approximation  for  X, 

(7)  X  ~7  • 

Hext  we  integrate  the  system  (3)  with  these  initial  conditions 

(8)  x(4)  =  -  1.80843 
u(4)  =  +  0.08 
x(4)  =  +  7-0 

on  the  interval  4<jts8  and  obtain  the  functions  XQ(t),  uQ(t),  \Q(t) 
on  that  interval. 

To  obtain  the  (n+l)  approximation,  after  having  calculated 
the  nth,  we  use  the  linearized  relations 

un+l 

-  0  un"  xn  *  (W  xn  X'  V.*.  '  l] 

*  (W  un  )  [-  <V  0]  +  (W-  »n  )  [-  (4  XX] 
0  , 


(9) 


n+l 


n+l 


'n+l 


together  with  the  three-point  boundary  conditions  of  Eqs .  (2). 
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The  results  of  a  numerical  experiment  are  summarized,  in  the 
Table  1. 


Table  1 

THE  FIRST  NUMERICAL  EXPERIMENT 


Initial 

Approx. 

Iteration 

Two 

Iteration 

Three 

Iteration 

Four 

True  Values 

x(k) 

-1.80843 

-1.80843 

-1.80843 

-1.80843 

-1.8084322 

u(4) 

+0.08 

+0.05644-54 

+0. 0794911 

+0.079366 

.079366909 

X(4) 

+7.0 

9.91541 

10.0004 

10.00000 

10.00000 

A  second  experiment  vas  carried  out  with  poorer  estimates  of 
the  initial  velocity  and  system  parameter.  The  results  are  presented 
in  Table  2. 


Table  2 

THE  SECOND  NUMERICAL  EXPERIMENT 


Initial 

Approx. 

Iteration 

Two 

Iteration 

Four 

Iteration 

Six 

True  Values 

x(4) 

-1.80843 

-1.8o843 

-1.80843 

-1.80843 

-1.8064322 

u(4) 

0.1000 

-2.0758 

0.599288 

.0792063 

.079366909 

A  (4 ) 

5.000 

3.87992 
_ 1 

11.4091 

9.99956  1 

...  . . .  .  .  i 

10.000 

A  third  trial  in  which  the  initial  estimate  of  the  system 
parameter  vas  taken  to  be  20  resulted  in  an  overflow,  so  that  no 
results  are  available. 

The  data  In  Bi.  (8)  were  generated  by  integrating  Van  der  Pol’s 


equation  with 
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X  =  10.0 

(10)  x(o)  =  1.0 

x(o)  =  0.0 

A  graph  is  presented  in  Fig.  1. 

The  FORTRAN  program  vhich  produced  these  results  is  given  next. 
The  subroutines  for  integration  of  a  system  of  ordinary  differential 
equations  and  the  numerical  solution  of  linear  algebraic  equations 
are 

Adams -Moulton,  Runge-Kutta  Integration  FAP  Coded  Subroutine 
TFQRTRAn),  Robert  Causey  and  Werner  L.  Frank,  Space  Technology 
Laboratories;  adapted  by  W.  L.  Sibley,  The  RAND  Corporation, 

June  1961,  RAND  7090  library  Routine  Number  X013.  (See  also 
SHARE  Distribution  #6<&.) 

Matrix  Inversion  vith  Accompanying  Solution  of  Linear 
Equations,  Burton  S.  Garbov,  Argonne  National  Laboratory, 
February  1959 >  70*+  SHARE  Routine  Number  664. 

Phase  II  in  the  program  is  used  vhen  the  number  of  observations 
equals  the  number  of  unknowns  and  phase  III  is  used  when  the  number 
of  observations  is  greater  than  the  number  of  unknowns. 

The  following  is  a  brief  summary  of  the  FORTRAN  code  for  phase 
II.  After  the  necessary  data  have  been  input,  the  initial  approxi¬ 
mation  is  generated  by  integrating  the  nonlinear  differential 
equations.  The  (n+l)st  approximation  is  obtained  by  integrating  the 
particular  and  homogeneous  equations,  determining  the  unknown  con¬ 
stants,  and  forming  the  linear  combination  of  the  particular  and 
homogeneous  solutions.  This  (n+l)st  approximation  is  printed  and 
stored  as  the  n'*'*1  approximation  in  preparation  for  the  next  iteration . 
If  the  initial  values  of  the  (n+l)st  approximation  agree  with  those 
of  the  n^*1  approximation  to  5  or  more  decimal  places,  the  problem  is 


considered  solved. 
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Solution  of  the  Von  Der  Pol  equation 
ir  X  =  10,  x(0)s0,  x(0)=1 
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CVDP021  VAN  DER  POL  -  PHASES  II.  Ill 

COMMON  SCRACH.T.Nl . NMAX » <MAX » HGR I D »NT IME » XOBS . W . H . P » A . B. X »U » Z » T I . 
1  PREV.N03S.IFLAG 

DIMENSION  SCRACH  ( 64 ) . T ( 1 50 ) . NT  I  ME ( 100 ) .XOBS (100) ,*<3,1001 ) . 

1  H(3»3»1G01) »P(3*1001 ) »A( 16.16) »8( 16.1) .PREV (3 ) 

C 

C  INPUT  NOBS  OBSERVATIONS  OF  X  AT  KNOWN  TIMES 

C  DETERMINE  SYSTEM  PARAMETER  LAMBDA  CZ  IN  FORTRAN) 

C  INPUT  AND  START 

1  CALL  INSTRT 
C 

C  K  ITERATIONS 

C 

DO  99  K= 1 »  KMAX 
C 

DO  2  1=1.150 

2  T (  I  ) =0 .0 
T ( 2 ) =T I 

T ( 3 ) =  HGR I D 
T  (  4  )  =  1  •  0 
T(8)*1.0 
T< 12 ) =1.0 
X=PREV( 1 ) 

U=PREV ( 2 ) 

Z=PREV(3) 

CALL  INT(T.12.N1.0«.0».0».0#.0«.0») 

N  =  1 
L  =  3 

DO  21  1=1.3 
DO  21  J= 1  ♦  3 
L  =  L  +  1 

21  H(I.J»N)=T(L) 

DO  22  1=1.3 

L  =  L  +  1 

22  P(I.N)=T(L) 

C 

C  INTEGRATE  P'S  AND  H'S 

DO  4  N=2 .NMAX 
X=W< 1 ,N) 

U=W ( 2  »N ) 

Z  =  W  (  3  » N  ) 

CALL  I NTM 
L  =  3 

DO  3  1=1.3 
DC  3  J=1 .3 
L  =  L+1 

3  H( I ,J,N)=T(L) 

DO  4  1=1,3 

L  =  L+1 

4  PC  I , N ) =T ( L ) 

C 

C  DETERMINE  CONSTANTS,  OR  INITIAL  VALUES 

CALL  CNSTNT 
T I ME=  T I 

PRINT  50,  K.TIME.CBd.l)  ,1  =  1.3) 

C 

C  NEW  VARIABLES 

DO  7  N=2 .NMAX 
DO  6  1=1,3 


n  n 
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W( I »N)=P( I ,N) 

DO  6  J= 1  *  3 

6  W( I »N)=W< I *N)+B( J»1)*H( J» I »N) 

FN=N-1 

T I ME=FN*HGR ID+T I 

7  PRINT  70,  TIME, (W(I,N) ,1=1*3) 

C 

C  COMPARE  CONSTANTS 

DO  8  1=1,3 

G=ABSF ( B ( 1,1) -PREVI I ) ) 

IFIG-.000001 ) 8,8,9 

8  CONTINUE 

GO  TO  1 

9  DO  10  1=1,3 

10  PREV (11=8(1,1) 

C 

99  CONTINUE 
GO  TO  1 
C 

50  FORMAT  ( 1H0//59X9H ITERATION,  1 3 //3.8X1HT  ,  1 3X 1HX  ,  1 9X1  HU ,  1 7X6HS YSTEM/ / 
1  30XF10.2,3E20.6) 

70  FORMAT (30XF10«2»3E2C,6) 

END 


CVDP022  INPUT-START,  VAN  DER  POL  PHASE  I  I  ,  I  1 1 

SUBROUTINE  INSTRT 

COMMON  SCRACH , T  »N1 *NMAX , KMAX ,HGR I D ,NT IME » XOBS , W , H , P , A , B , X , U , Z , T I , 
1  PREV, NOBS, IFLAG 

DIMENSION  SCRACH (64) , T ( 1 50 ) , NT  I ME < 100 ) ,XOBS ( 1 00 ) ,W(3,1001 ) , 

1  H(3»3*1001),P(3»1001),A(16»16),B(16,1> ,PREV<3) 

C  INPUT 

READ  110, Nl» NOBS, KMAX ,NMAX , HGR I D , (NT  I ME (I ) ,XOBS( I ) ,I  =  l,NOBS) 

I F ( NMAX ) 9 , 9 , 1 

1  PRINT  10, N1 , NOBS, KMAX, NMAX, HGRID, (NTIME( I ) »XOBS( I ) ,1=1, NOBS) 

START 

iflag=i 

READ  120 »TI »X»U»Z 
DO  2  1=1,150 

2  T( I ) =0,0 
T ( 2 ) =T I 
T ( 3 ) =HGR  I D 
T ( 4 ) =X 
T ( 5 ) =U 
T  ( 6  )  =2 

CALL  INT(T»3»N1,0.»0.»0.»0.,0.*0.) 

K  =  0 

PRINT  50,  K,T(2),(T(I) ,1=4,6) 

C 

DO  4  N=2 ,NMAX 
CALL  INTM 
DO  3  1=1,3 

3  W( I ,N)*T( 1+3) 

4  PRINT  70,  T(2 ) , (T( I ) , 1=4,6) 

C 

I FLAG=2 
PREV ( 1 ) =X 
PREV(2)=U 
PREV(3)=Z 


c 


RETURN  ** 

9  CALL  EXIT 
10  FORMAT (1H130X 

1  58HVAN  DER  POL  -  PHASE  II  -  DETERMINATION  OF  SYSTEM  PARAMETER// 

2  20X4 1 10, Flu. 4/ < 1 8X 1 1 2 * El 6 .6 » 1 12 »E16. 6. 1 12 »E16 .6 . 1 12 »E16 . 6 ) ) 

50  FORMAT (1H0//59X9H ITERATION. 1 3 / /38X1HT . 13X 1HX » 19X1HU » 1 7X6HSYSTEM// 
1  30XF10.2.3E20.6) 

70  FORMAT  (30XF10.2.3E20.6) 

110  FORMAT (4I5*F10«2»2( I4.E11.6 )/ ( 4  ( I4.E11.6) ) ) 

120  FORMAT (4E12. 6 ) 

END 


CVDP023  DAUX  -  VAN  DER  POL  -  PHASE  1 1 • 1 1 1 

SUBROUTINE  DAUX 

COMMON  SCRACH  » T  »N1  .NMAX.KMAX  »HGR  ID  »NT  IME  »XOBS  .vJ.H.P.A.B.X.U.Z.TI 
1  PREV»NOBS» IFLAG 

DIMENSION  SCRACH ( 64) *  T  « 150 ) » NT  I ME ( 100 ) ♦ XOBS ( 1 00 ) *W(3»1001 ) * 

1  H( 3.3.10C1 ) *P<  3*1001 ) »A( 16. 16) »B( 16»1 ) »PREV(3) 

C 

GO  TO  (1*2)  .IFLAG 

C  NONLINEAR  EQUATIONS 

1  T ( 7 ) =T ( 5 ) 

T(8)s-(T(4)**2-1. )*T(5)*T(6)-T(4) 

T ( 9 ) =0  « 

RETURN 

C  LINEAR  EQUATIONS 

2  XX=X*#2 
AA=-2.*X*U*Z-1. 

BB=Z*(1.-XX) 

CC=U*(1.-XX) 

L  =  1 3 

DO  3  1=1.4 
L  =  L+3 

T ( L ) =T ( L-l 1 ) 

T(L+1 )=AA*T(L-12 )+BB*T(L-ll )+CC*T (L-1C) 

3  T ( L+2 ) =0  •  0 

T (L+l ) =T ( L+l ) +U*Z*( 3.*XX-1. ) 

RETURN 

END 


CVDP024  CNSTNT  -  PHASE  II 

SUBROUTINE  CNSTNT 

COMMON  SCRACH,  T.Nl  .N.MAX  ,  KMAX  .  HGR  I D  .NT  I  ME  .  XOBS  »  W  »H  .  P  .  A  .B .  X  ,  U  .  Z  » T  I 
1  PREV  »N03S  » IFLAG 

DIMENSION  SCRACH (64) , T ( 1 50 ) » NT  I ME ( 1 00 ) .XOBS ( 100) »W(3»1001) » 

1  H(3»3,10Gl),P(3.1w01),A(16.16),o(16,l)  .PREV (3) 

C 

DO  1  1=1,3 
N  =  NT IME (  I  ) 

B ( I ,1 )=XOBS( I )-P( 1 ,N) 

DC  1  J=1  ,3 

1  A( I ,J)=H< J.l.N) 

CALL  MATINV(A,3»B.1»DET) 

DC  2  1=1,3 

2  W( I ,1 )=B( I .1) 

RETURN 

END 
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VII.  COMPUTATIONAL  CONSIDERATION  -  II 

The  successive  approximation  method  proposed,  which  is  abstractly 
equivalent  to  Nevton's  method  for  extracting  roots,  shares  with  that 
procedure  the  desirable  property  of  quadratic  convergence.  This 
means,  roughly,  that  the  number  of  correct  digits  is  doubled  with 
each  additional  iteration.  Furthermore,  the  computational  load  at 
each  stage  is  light,  requiring  only  the  integration  of  some  initial- 
value  problems  and  the  solution  of  some  linear  algebraic  equations . 
Nevertheless,  difficulties  can  arise.  Let  us  discuss  some  of  these. 

In  the  first  place,  a  solution  to  a  nonlinear  multipoint  boundary- 
value  problem  need  not  exist,  nor  need  it  be  unique.  With  well- 
formulated  problems  arising  from  physical  sources,  however,  we  would 
not  expect  this  to  be  a  source  of  difficulty.  Of  more  practical 
import  is  the  fact  that  if  the  initial  approximation  is  too  far  re¬ 
moved  from  the  solution,  the  iterations  may  not  converge.  This 
possibility  deserves  further  investigation. 

We  have  found  that  the  integration  of  the  initial  value  problems 
can  be  difficult,  for  we  wish  to  choose  a  grid  size  that  is  suffi¬ 
ciently  small  to  give  the  required  accuracy  yet  not  so  small  as  to 
involve  excessive  computing  costs.  Furthermore  the  numerical  solu¬ 
tion  of  the  linear  algebraic  equations  for  the  multipliers  of  the 
solution  of  the  homogeneous  equations  can  be  difficult,  for  though 
the  matrix  involved  is  nonsingular,  it  can  be  ill-conditioned.  A 
promising  method  for  overcoming  this  difficulty  is  described  in  [13]. 

The  successive  approximation  scheme  proposed  involves  storing 
the  values  of  the  dependent  variables  at  one  stage  to  calculate  their 


1 6 


values  in  the  next.  The  high-speed  memory  limitations  of  the  com¬ 
puter  being  used  can  he  exceeded.  A  method  for  overcoming  this 
is  described  in  fl4]. 
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